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Cramer-Rao Lower Bound for Point Based Image 
Registration with Heteroscedastic Error Model for 
Application in Single Molecule Microscopy 
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Abstract —The Cramer-Rao lower bound for the estimation 
of the affine transformation parameters in a multivariate het¬ 
eroscedastic errors-in-variables model is derived. The model is 
suitable for feature-based image registration in which both sets of 
control points are localized with errors whose covariance matrices 
vary from point to point. With focus given to the registration of 
fluorescence microscopy images, the Cramer-Rao lower bound for 
the estimation of a feature’s position (e.g. of a single molecule) 
in a registered image is also derived. In the particular case 
where all covariance matrices for the localization errors are 
scalar multiples of a common positive definite matrix (e.g. the 
identity matrix), as can be assumed in fluorescence microscopy, 
then simplified expressions for the Cramer-Rao lower bound are 
given. Under certain simplifying assumptions these expressions 
are shown to match asymptotic distributions for a previously 
presented set of estimators. Theoretical results are verified with 
simulations and experimental data. 

Index Terms —Image registration, Cramer-Rao lower bound, 
generalized least squares, fluorescence microscopy. 


I. Introduction 

I MAGE registration is the process of overlaying two or more 
images of the same scene m . Image registration techniques 
can be divided into two categories; intensity-based registration 
where gray scale values are correlated between images, e.g. 
El a, 0, and feature-based registration, whereby corre¬ 
spondence between the two images is determined through the 
matching of distinct features common in both images e.g. Q, 
©. 

This project is motivated by an important problem in 
single molecule microscopy, a recent major advancement in 
fluorescence microscopy which allows individual fluorescently 
labeled molecules to be imaged using optical microscopy 
techniques and individually localized with accuracies in the 
very low nanometer range El, El, IS. In a typical experiment 
two different proteins in a cell are labeled with different 
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Fig. 1: A diagram illustrating fiducial markers visible in 
both images (left). Each fiducial marker’s position is located 
with a measurement error. Registration requires finding the 
transformation that best aligns the fiducial markers (right) with 
respect to an appropriate minimization problem. 


fluorescent markers. The biological information is obtained 
from the relationship between the two labeled sets of proteins. 
The imaging experiment consists of taking one exposure 
for each of the labeled proteins, often using two cameras, 
each equipped with a wavelength dependent optical filter to 
capture the emission of the fluorescence for the corresponding 
proteins. In this fashion we obtain two different images each 
displaying different aspects of the sample. In order to analyze 
these images they need to be registered, as it cannot be 
assumed that the cameras are aligned to the degree that is 
necessary to guarantee the nanometer level accuracy which is 
required to obtain the appropriate information. Registration is 
typically achieved by incorporating fiducial markers, usually 
small nanometer size beads, into the sample whose fluorescent 
properties are such they can be imaged in both cameras. 
These fiducial markers can therefore serve as control points 
(CPs) for feature-based registration. The characterization of 
the registration errors is critical in assessing the deterioration 
of the localization accuracy of a single molecule due to the 
registration. A number of further single molecule microscopy 
experiments lead to the same underlying registration problem. 
One important such example arises from the correction of drift 
in time lapse experiments. 
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Previous statistical studies on CP registration Qoi, imiii, 
ns, m, insi, m, ini assume one set of positions for the 
control points is taken as truth and errors exist in only the 
second set of positions. In such examples a multivariate linear 
regression is used for the data model and the set of localization 
errors are sometimes referred to as the fiducial localization 
errors (FLEs) iflTl . ifT^ . In contrast, a recent study ifTSl pre¬ 
sented a multivariate errors-in-variables (EIV) formulation of 
the control point image registration problem important to the 
microscopy application. The errors-in-variables formulation is 
necessary to model the situation when a ground truth for the 
CP locations is unavailable and the errors in measuring the CP 
locations are present in both images that are to be registered. 

Importantly, the microscopy application dictates that a het- 
eroscedastic model is used as the measurement errors have to 
be assumed to have different covariance matrices for different 
CPs. A central aspect of the registration problem is the 
estimation of the registration transformation. There have been 
attempts in im, lEO) to estimate registration parameters for 
heteroscedastic errors under an EIV model when the transfor¬ 
mation is assumed rigid (rotation and translation only) with the 
heteroscedastic EIV (HEIV) algorithm; an iterative procedure 
that finds an optimal solution to the HEIV model. However, 
estimator distributions were only determined through boot¬ 
strapping methods. In the microscopy setting we take the 
more general assumption that the registration transformation 
be affine, allowed due to the high geometrical precision of 
modern microscope objectives. In ifTSl we have shown that for 
this data model a generalized maximum likelihood estimator 
is equivalent to a generalized least squares estimator. Using 
prior results we were able to obtain asymptotic results on 
the distributions for the estimators for the transformation 
parameters. Eor a specific heteroscedastic noise model where 
covariance matrices are scalar multiples of a known positive 
definite matrix, closed form expressions for estimators of the 
affine transformation parameters were derived. This particular 
model is applicable in a fluorescence microscopy setting where 
fiducial markers (e.g. fluorescent beads) act as the CPs but are 
each localized with differing degrees of accuracy. 

Registration performance is typically quantified by the 
fiducial registration error (ERE), which is the root mean- 
square distance between fiducial markers after registration, and 
most importantly the target registration error (TRE) which is 
the difference between corresponding points (other than the 
fiducial markers) after registration. The distribution of the TRE 
has been of much interest. Under the multivariate linear regres¬ 
sion model (which as stated is inappropriate in fluorescence 
microscopy) Eoi, im derive approximate expressions for the 
root mean square of the TRE’s absolute value and ifT^ gives 
its approximate distribution in the case where the registration 
transformation is assumed rigid (rotation and translation only) 
and ELEs are independent and identically distributed (iid) 
zero-mean Gaussian. Readers interested in the effect of biased 
ELEs are directed to inj. Anisotropic iid FLEs are first 
considered in na and M derives the maximum likelihood 
estimators for the rigid transformation parameters along with 
the associated Cramer-Rao Lower Bounds on their variance for 
this model. Heteroscedastic FLEs are considered in OH and 


using a spatial stiffness model they derive expressions for the 
root mean square TRE. An overview of these methods is given 
in lEa, together with procedures for the optimal selection of 
fiducial markers with respect to minimizing the TRE. 

In 1181 asymptotic distributions were found for the TRE 
under the multivariate errors-in-variables data model and 
affine transformation assumption required for fluorescence 
microscopy. Eurther to this, in lITSl the asymptotic distribution 
was also found for the localization registration error (ERE), 
a newly defined measure of registration error that combines 
both a localization error and the TRE of a feature (e.g. single 
molecule) that is not used in the registration. 

The quality of a single molecule experiment is assessed by 
the accuracy with which single molecules are localized in the 
particular experiment 12^ . Here the localization accuracy is 
interpreted as the standard deviation of an unbiased location 
estimator ll24l . In Il24l and ll25l the fundamental limit of 
localization accuracy was introduced as the Cramer-Rao lower 
bound (CRLB) for the location estimation problem, in the 
context of ideal experimental conditions such as an infinite size 
photon detector without pixilation artefacts and without other 
extraneous noise sources. This measure has proved a reliable 
predictor for the best possible accuracy that can be achieved 
with a specific single molecule experiment ED, Ez). 

Due to the importance of registration in single molecule 
experiments the question therefore arises how the uncertainty 
introduced during the registration process influences the local¬ 
ization accuracy for a single molecule that has been registered. 
To this end, a major aspect of this manuscript consists of the 
derivation of the CRLB for the registration problem for several 
data models that are of relevance here. 

The CRLB has been derived for registration problems 
before. The work of ll28l and E9l consider the CRLB for 
feature-based and intensity-based registration performance in 
several scenarios of more general affine transformations be¬ 
tween the two images, as well as a polynomial based non¬ 
linear transformation. However, they restrict themselves to the 
homoscedastic case, i.e. when all CP measurement errors have 
equal covariance matrix and consider only the CRLB of the 
transformation parameters themselves. 

This paper provides the CRLB for registration performance 
when a general affine transformation is assumed and in the 
case of heteroscedastic CP measurement errors assumed zero- 
mean and Gaussian. We give particular focus to a fluorescence 
microscopy setting and not only consider the CRLB in esti¬ 
mating the transformation parameters, but place emphasis on 
finding the lower bound for the covariance matrix of the LRE, 
a concise and informative measure of registration performance. 
The square root of the diagonals of this covariance matrix 
(the standard deviation of the LRE in each dimension) is 
the accuracy with which single molecules are localized post¬ 
registration in each dimension. 

In Section HI] we formulate the registration problem and 
define the LRE as introduced in IlSl . In Section uni we 
derive the CRLB for the affine transformation parameters in 
the most general heteroscedastic setting. In Section |IV| we 
derive the CRLB for estimating the unknown position of a 
feature in the registered image, in turn giving a lower bound 
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on the variance of the LRE. In Section [V] we consider the 
specihc case when the covariance matrices of the measurement 
errors are a scalar multiple of a common matrix as is the 
case in a fluorescence microscopy setting. In fluorescence 
microscopy it is reasonable to approximate the covariance 
matrices as multiples of the identity matrix. Further, it is 
common for the affine transformation matrix to be a scalar 
multiple of a unitary matrix (a combination of scaling, rotation 
and reflection). In such a scenario we derive an explicit 
expression for the lower bound of the covariance matrix of 
the LRE that reveals a more intuitive view of these complex 
expressions. Importantly, this expression is identical to that 
of the asymptotic covariance of the LRE. This was derived 
in GSl under equivalent assumptions when the registration 
parameters are estimated using the corresponding generalized 
least-squares or equivalent maximum likelihood estimator. In 
Section |Vl] we verify the theory with simulation studies and 
show the generalized least-squares estimator of m attains 
this lower bound. We conclude by considering real microscopy 
imaging data and show the CRLB results presented in this 
paper are appropriate in an experimental setting. 

II. Formulation 

We consider the registration experiment formulated in ifTSl . 
There are K CPs located in both image 1, denoted Ii C 
and in image 2, denoted I 2 C (d = 2 or 3). These 
CPs have true locations {xi^k G Ii,k = and 

{x 2 ,k G 2^2, fc = respectively, and these CP coor¬ 

dinates are related by the affine transformation T : 
where X 2 ,k = T{xi^k) = Axi^k s, k = l,...,iT, with 
invertible A € and s G R"^. The true positions of 

the CPs are not known in either image and instead must be 
measured with errors. We therefore observe the CP locations 
as { 1 / 1 ,fc €li,k = and {1/2,fe G l 2 ,k = 

where yj^k = Xj^k + ^j,k, k = j = 1,2. The 

term ej^k G R'^ is a random measurement error, sometimes 
referred to as the fiducial localization error (FEE), and are 
each assumed zero mean and to have individual covariance 
matrix 17^ ^ > 0 (where we use notation M > 0 if matrix M is 
positive definite and M > 0 if it is non-negative definite). All 
measurement errors are assumed to be pairwise independent 
across the CPs. 

Let us define the R'^^^ matrices Xj = [xj^i, ...jXj^k], 
Yj = Sj = j = 1)2. 

The measured control point locations can be conveniently 
represented as Yi = Xi+ £i and y 2 = X 2 + £ 2 - The latter of 
these can be equivalently represented as Y 2 = AXi-\-s\^+£ 2 , 
where ^ is the matrix transpose and 1^ is a column vector 
of length K with every element taking the value 1. If we 
further define the stacked matrices X = 

Y = \Yi ,Y 2 \'" and £ = then the system of 

equations can be condensed into the single matrix equation 

y = AXi-f al^-f £•, (1) 

where a = [0^, and A = [Id, with Id representing 

the d-dimensional identity matrix. The columns of £ are 


independent with fcth column Ck = [^i kX 2 kl'^ having mean 
zero and known positive definite covariance matrix 


rife = covjefc} 


Hi.fe 0 

0 172.fe 


( 2 ) 


where cov{u} denotes the covariance matrix of a random 
vector V. 

Models of type O are called errors-in-variables models. 
When covariance matrices {£tk,k = l,...,Ar} all equal the 
same matrix 17o > 0 we have a homoscedastic errors-in- 
variables model. Under the homoscedastic assumption O 
is equivalent to the registration formulation of the CRLB 
study by 1281 and 1291 . When the K covariance matrices 
{£lk,k = 1, ...,K} are in general not equal then we have a het- 
eroscedastic errors-in-variables model. It is the heteroscedastic 
assumption that this study focuses on. 

Image registration requires estimating the transformation 
parameters A and s whose elements we can represent in the 
transformation parameter vector 6t = [vec(A)^, s^]^. In the 
strict homoscedastic case l30l defines the generalized least 
squares (GLS) estimators of A and s and shows them to be 
equivalent to the maximum likelihood (ML) estimators under 
the assumption of CP measurement errors being Gaussian. 
Further to this, closed form expressions for the estimators 
of A and s are given along with their joint and marginal 
asymptotic distributions. The work of m considers the most 
general heteroscedastic model, where under the assumption 
of Gaussian measurement errors the maximum likelihood 
estimators for A and s are presented along with an iterative 
method for their computation and their joint and marginal 
asymptotic distributions. Recently in ESI, a heteroscedastic 
generalized least squares estimator is defined in an extension 
to the homoscedastic formulation of 1^ and is shown to be 
equivalent to the maximum likelihood estimator considered 
in ED. In the special case where covariance matrices for 
the measurement errors are of the form flk = Vk^o, where 
rik G R"*" and Ug > 0 — termed the weighted covariance 
model — then M gives closed form expressions for the 
estimators of A and s and determines their joint and marginal 
asymptotic distributions. This in turn is used to give concise 
expressions for the first and second moment of the TRE and 
LRE, measures of registration error that we now formally 
define. 


A. Registration Errors 

As has been stated in Section H] the TRE is a commonly 
used measure of registration performance. Here we give its 
definition when the registration transformation T is assumed 
affine with matrix parameter A and vector parameter s (see 
Figure |2]i. 

Definition II.l. Let A and s be the registration transformation 
parameters and let A and s be their respective estimators. The 
target registration error (TRE) r : Ii —>■ for an arbitrary 
point xi G Xi with corresponding mapped position in I 2 of 
X 2 = Axi -\- s is defined as t(xi) = X2 — (Axi -f s) = 
Axi -I- s — {Axi + s). 
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Image 1 Image 2 

-, I-, 



I >1 • 

L_' L_' 

Fig. 2; A diagram illustrating the target registration error as 
the difference between the true position of an arbitrary point 
in I 2 (black/upper cross) and its registered position using the 
estimated registration parameters (green/lower cross). 


Image 1 Image 2 



Fig. 3: A diagram illustrating the localization registration 
error as the difference between the true position of feature 
in I 2 (black/upper cross) and its registered position using its 
estimated position in Ii (green/lower cross) and the estimated 
registration parameters. 


The LRE is defined in |[T3 and is of particular use in 
fluorescence microscopy registration experiments. Suppose we 
have a feature (e.g. a single molecule) that is visible in image 
1 but not in image 2 (and therefore is not involved in the 
registration process). The LRE gives the error with which it 
is localized in image 2 after registration (see Figure O. 

Definition II.2. Let A and s be the registration transformation 
parameters and let A and s be their respective estimators. 
For a feature (e.g. single molecule) in Ii with true and 
measured locations xi^p ond yi^F = respectively, 

the localization registration error (LRE) If is defined as the 
difference between the true position in I 2 . given by X 2 ,f = 
Axi,f + s, and the registered position X 2 ,f = + s, i.e. 

If = X 2 ,F — X 2 ,F- 

The standard deviation of an element of £f is the accuracy 
with which a feature/single molecule can be localized in 
that dimension post-registration. The covariance matrix for 
£f, denoted is identical to the covariance matrix for 
the estimator X 2 ,f and therefore the CRLB for estimating 
X 2 ,f is a lower bound for fig, i.e. if we denote the CRLB 
matrix for estimating X 2 ,f as Cff then fig > Cff (where 
notation Lli > Cff means flu — Cff > 0, i.e. fli — Cff is 
non-negative definite). Further discussion on the relationship 
between the TRE and LRE can be found in M- 

III. CRLB FOR AFFINE TRANSFORMATION PARAMETERS 

The vector of unknown parameters is given as 9tc = 
where 0t = [vec(A)^, is the +‘^) 

vector of affine transformation parameters and 9c = 


\x\ kV i^ *^he vector of CP location pa¬ 
rameters. We make the assumption that measurement errors 
are independent, are zero-mean (validated in ll26l l and are 
Gaussian (validated in Section FVll i with covariance matrix flk 
of form (|2]l, i.e. ^ A^2d(0, flfc). The likelihood function is 
therefore given as lfT8l . ll3T1 

L{9Tc\yi, yK) = (^11 

X exp - yk)'^flf^(yk - , 

where = [x"[f.,x'^ fy’", |f2fc| denotes the determinant of 
flk and yi,...,yK are the stacked vectors of measured CP 
locations. The corresponding log-likelihood is 


Ci9Tc\yi, ■■■, yK) = 

Kd 1 ^ I ^ 

^ lii( 27 r) —— In —/Tfc) {]Jk~yk)- 

k=l k=l 


It is well established that for the multivariate normal distri¬ 
bution with covariance that is independent of parameters the 
Fisher information matrix (FIM) for the parameter vector 9tc, 
denoted J{9tc), is given as ll^ p. 47] 


J{Stc) = 


^ dyl duk 


(3) 


and the CRLB matrix is given as C{9tc) = J With 

the CRLB matrix denoted as 


C[9tc) 


Ctt Ctc 
C cT Ccc 


(4) 


the diagonals of Ctt are the CRLBs for estimating the 
transformation parameters and the diagonals of Ccc are the 
CRLBs for estimating the control point locations. We are 
therefore primarily interested in the diagonals of Ctt- P is 
shown in Appendix IaI that 


K 

J{9tc) = X! 


HlflflHk HlflflCk 


where Fk = 0 h, Gfe = ® A, Hk = 

Id C' X-^ Id 


ik}\T 


Here, Cp represents the jth standard basis 
vector of R^, (i.e. vector of length p with 1 placed in the jth 
entry and zeros everywhere else) and (g> denotes the Kronecker 
product. It follows from the block matrix inversion of J{9tc) 
that 

-1 


Ctt = [Shh — Shg (Sff + Sac) ^ Sj)c^ 


(5) 


where Shh = Eti Hl^^kHk, Shg = Ef=i Hl^flGk, 
Sff = Sfc=i Fkfli kFk and Sgg = G]^fl 2 ,kGk- This 
is the CRLB for the transformation parameters in the most 
general heteroscedastic errors-in-variables model considered 
in on . Equations dm . (fT3]) and (fT4l) in Appendix IaI pro- 
vide expressions for the sums Shh, Shg and Sff + Sgg, 
respectively. 
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IV. Feature localization 

Let us now consider including the localization of a feature 
(e.g. single molecule) into the expression. For this we include 
the term yi,F - the observed location of the feature in Xi, and 
the unknown parameter X 2 ,f - the true position of the feature 
in I 2 that we wish to estimate. The associated localization 
error has covariance matrix As previously stated in 

Section III-Al the LRE If G is the difference between 
the estimator X 2 .f = Mh-X + s and the true value X 2 ,f and 
hence the CRLB for estimating X 2 ,f provides a lower bound 
for the covariance matrix of £f- 

The combined log-likelihood function for all parameters 
Oftc = where 9f = X2,f^ given the observed 

data is now 

CiOFTclvi, ■■■,yK,yix) = - 2 

- 2 “ ^'ixV^ixiy^X - Mi,f)- Y 

1 iC ^ K 

k=l k=l 

In terms of the unknown parameter Of = X 2 ^f we can write 
Pi F = A“^(a; 2 ,F ~ s). The FIM for the complete parameter 
vector Oftc is shown in Appendix iBl to be given as 


V. CRLB EXPRESSIONS FOR WEIGHTED COVARIANCE 
MODEL 

Section ITVi provides the CRLB for localizing a feature/single 
molecule in the most general heteroscedastic registration 
model. While these results provide a very general solution 
to our problem, we will now investigate special cases that 
are of interest in their own right through their relevance 
in applications. In addition, in these special cases we can 
obtain significant simplifications of the above expressions that 
provide useful insights for experimental design considerations. 
In this section we look to the weighted covariance model 
formulated in Ea, in which we make the assumption that 
covariance matrices for the measurement errors are of the form 
= Vk^o where rjk S K’*' and flo > 0, for all k = 1, ..., K. 
Here, we consider the following further assumption. 
Assumption I. Covariance matrices have the forms: Hi.o = 
(Ji qI 2 , ^ 2,0 = o ’2 0 ^ 2 . = Cl pl 2 , and transforma¬ 

tion matrix A = where R is a unitary matrix (rota- 
tion/refiection) and € R."'' is a scaling factor. 

The transformation vector s is arbitrary. The assumption here 
that the covariance matrices are some scalar multiple of the 
identity matrix is a reasonable assumption in fluorescence mi¬ 
croscopy and exact in the case of a non-pixelated detector ll24ll . 
The assumption on the transformation matrix is a common 
type of transform experienced in registration. 


J{0ftc) = 

UpT 

0 

where Dtt = and Dpr = D'f^iSp^t, with 

A~^ [xl^F ® Id, Id] ■ 

Representing the inverse FIM of Oftc as 

CpF CfT CfC 

Ctf Ctt Ctc , 

CcF CcT Ccc 

it is shown in Appendix that the sub-block [ ] is 

identical to C{0tc) in @ (as one would expect from the 
fact that the feature/single molecule is not involved in the 
registration process), and the CRLB matrix for estimating 
X2,f, the location of a feature/single molecule in the registered 
image, is given by 

Cff = {A~'^Vtf\,A~^ — Dft{Dtt + ) 

( 8 ) 

where Ctt is the CRLB matrix for estimating the transforma¬ 
tion parameters given in (|5]l. The d diagonal elements of Cff 
are the CRLBs for estimating the respective elements of X2,f, 
and with Vli > Cff offers the lower bounds on the variances 
of the respective elements of the LRE tp- 


C{0ftc) = J ^{Oftc) = 


Dp = 


Dt = 


dp 


1,F 


do]; 


dp 


1,F 


dOl 


Dft 

Dtt + Shh 




HC 


0 

Shg 

Sff+Sgc 


, (7) 


A. CRLB for estimating transformation parameters 

Let us define the following quantities that will be used here; 
7 = (l/Ar)Ef=i»7fe\ Xk = (k = E = 

i^/K)J2k=iVk^Xk, Xt,k = 62 *^ ®xl,. (* = 1,2 and k = 
xi = {l/K)Yj^^i'q'^^Xi.k, X, = 0xf (i = 
1 , 2 ), = 5 - j-^xixf and F, = (i = 1 , 2 ). 

Under Assumption I it is shown in Appendix |D] that the 
expression for Ctt in © simplifies to 


Ctt = 


2 ,..2 
K 


^2,0; 


X 


0 -Tj 

0 -rj_ 

-Fi -Fs j-^h+l-^TiXl + T2XI) 


(9) 


B. CRLB for estimating the location of a feature/single 
molecule in the registered image 

1) General model: Under Assumption I it is shown in 
Appendix |E] that the CRLB matrix for estimating the location 
X 2 x of a feature/single molecule is given as 



xIf 1 






Xk 

0 


\ 

0 

Xk 

Xlk 


. Xi,k 

X2,k 

-^2 

J 


X [xJf,^'^ ■ (10) 
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2) Simplified model: Let us consider the case where rjk is 
independent of CP position and CPs are centrally and symmet¬ 
rically distributed in the image space. This model leads to the 
following set of assumptions that are appropriate for large K 
and asymptotically exact. These assumptions naturally arise, 
for example, when considering the experimental disposition of 
fluorescent beads in a specific microscopy experiment setting, 
where the beads can be assumed to have a circular Gaussian 
spatial distribution iTT^ . 

Assumption II. Approximate {1/K) 'Yl!k=i Vk^Xk = k'^l 2 kind 

iyK)j:t,Vk^X,,k = 0 . 

Under Assumption I and II it is shown in Appendix |F] that 
the CRLB matrix for estimating the feature/single molecule 
location X 2 ,f is given as 



3) Applying to a microscopy setting: Consider a fluores¬ 
cence microscopy example were Assumptions I and II are 
satisfied. That is, we have the weighted covariance model with 
rjk = Xfil, where Ni^k is the number of photons associated 
with control point k in Ii, k = and 

where o-f q = Ci cr| g = C 2 /C 133 . Constant Q, j = 1, 2, 
is a known localization accuracy parameter associated with Xj 
and is a function of the numerical aperture, photon wavelength 
and point spread function (see lUS] p. 6296]). Constant c is 
the constant of proportionality assumed in ifTSl to exist such 
that N 2 ,k = cNi^k where N 2 ,k are the number of photons 
associated with the kth CP in l 2 - This gives 


Xlk — 


1 

^k 


Cih 0 \ 

0 C-1C2/2 ;■ 


In this situation we have 7 = iVi where Ni is the mean 
photon count for the CPs in Xi. Therefore the CRLB matrix 
for estimating X 2 ^f 7 the location of the single molecule in X 2 , 
is given as 


CfF — kJi f^2 


1 




Photon counts are independent of CP position and there¬ 
fore under the (asymptotically exact m) assumption that 
= Nik^, where (l/AT) Xfe = (a measure of 

the spread of the CPs), then 

u, > Cff = ealFh + (^'1^ + ^) (1 + ^2- 

( 11 ) 

This expression for the CRLB in estimating X 2 ,f, and hence 
the lower bound for exactly matches the large K ex¬ 
pression for found in ifTSi p. 6297] when the generalized 
least squares estimator for the weighted covariance model is 
used. The fundamental limit (i.e. the theoretical lower bound) 
of localization accuracy for a single molecule in a pair of 
registered images is therefore bound by a term that depends 
on K (the number of CPs used in the registration process) 
and their associated photon counts, along with that gives a 
measure of the spread of the CPs in the image. We note there 
is no dependence on the translation parameter s. 


VI. Simulations and experimental verification 

In this section we verify the CRLB results given in this 
paper with computational simulations and a real data experi¬ 
ment. 

A. Simulation Studies 

In these computational simulation studies we consider a 
microscopy experiment where we register a pair of different 
coloured monochromatic images each captured using an op¬ 
tical system with identical numerical aperture and point-spread 
function. The measurement error in localizing the fcth CP in 
jth image Xj (j = 1, 2 ) has zero mean and covariance matrix 
iC]/Xj,k)l 2 where Q = A|^^/(dTr^nl.) The photon 

wavelength Xj^em associated with each image is 540nm and 
650nm respectively, uf is, the numerical aperture and assigned 
a typical value of 1.4 and Nj^k is the photon count associated 
with the /cth control point in the jth image. 

1) Rotation: CPs are arranged in a square grid of side 
length 81/im in the object space with varying numbers of 
points within that grid, and therefore K is restricted to 
the square numbers from 4 to 64. The photon counts as¬ 
sociated with each control point are observed realisations 
of a uniformly distributed random integer on the inter¬ 
val [5000,10000]. In Xi is a single molecule at position 
(16/rm,20/rm) from the center, with which a photon count of 
1000 is associated. Affine transformation matrix A is a rotation 
matrix of angle 30 degrees and affine transformation vector s 
is [4.8/rm,4.8^m]^. 

We look to verify the CRLB for estimating the trans¬ 
formation parameters as given in dUl and the CRLB for 
estimating the position of the single molecule in the registered 
image (equivalently the lower bound of Vli) in (ITTl l. This is 
achieved by estimating the transformation parameters using the 
generalized least squares estimator for the weighted covariance 
model as developed in lITSl . The empirical standard deviations 
of interest are computed using 10 ® simulations and shown in 
Figure |4] 

2) Shear: CPs are arranged in a square grid of side length 
81/im in the object space with iC = 9. The photon counts 
associated with each control point are observed realisations 
of a uniformly distributed random integer on the interval 
[5000,10000] and covariance matrices for the measurement 
errors are of the same form as in Section VI-Al. In Xi is 
a single molecule at position (16/im,20/im) from the center, 
with which a photon count of 1000 is associated. Affine 
transformation matrix A is a shear matrix of type A = (J ^) 
where shear parameter A is varied between values of 0.1 and 
0.9. Transformation vector s is [4.8/im,4.8/im]^. 

We look to verify the CRLB for estimating the trans¬ 
formation parameters as given in (Ell and the CRLB for- 
estimating the position of the single molecule in the registered 
image (equivalently the lower bound of ili) in (ITTll . This is 
achieved by estimating the transformation parameters using the 
generalized least squares estimator for the weighted covariance 
model as developed in lITSl . The empirical standard deviations 
of interest are computed using 10 ® simulations and shown in 
Figure |5] 
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Shear Parameter A. 
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0.95 
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Shear Parameter A 


Fig. 4; The line indicates the square root of the CRLB and 
the crosses mark the sample standard deviation plotted as a 
function of the number of CPs K for the following estimated 
parameters: (a) the first dimension of the unknown single 
molecule position parameter X 2 ,f, (b) si. the first element 
of transformation vector s, (c) an, the leading element of 
transformation matrix A, (d) 021 , the (2, l)th element of A. 
In this simulation study CPs are in a grid conhguration and 
the transformation consists of a rotation and translation (see 
Section IVI-Al I for more details). The vertical axes in (a) and 
(b) have units as nanometers, the vertical axes in (c) and (d) 
are unitless and are on the scale xl0“^. Results are based on 
a sample of 10® simulations. 


Fig. 5: The line indicates the square root of the CRLB and 
the crosses mark the sample standard deviation plotted as a 
function of shear parameter A for the following estimated 
parameters: (a) the hrst dimension of the unknown single 
molecule position parameter X 2 ,f, (b) si, the hrst element 
of transformation vector s, (c) an, the leading element of 
transformation matrix A, (d) 021 , the (2, l)th element of A. In 
this simulation study there are nine CPs in a grid conhguration 
(see Section rVl-A2l for more details). The vertical axes in (a) 
and (b) have units as nanometers, the vertical axes in (c) and 
(d) are unitless and are on the scale x 10“^. Results are based 
on a sample of 10® simulations. 


3) Asymptotic covariance versus CRLB: It has been men¬ 
tioned in Section IV-B3I that under Assumption I and II the 
lower bound for in ([n]l matches the large K covari¬ 
ance matrix of the LRE given in m when transformation 
parameters are estimated using the generalized least squares 
estimator. We now consider relaxing Assumption I such that 
rio is no longer the identity matrix and look at how the CRLB 
for estimating X 2 ,f compares with the more general large K 
covariance matrix expression in ifTSl p. 6295]. 

We have exactly the same experimental set-up as in Section 
IVI-All except the measurement error in localizing the fcth 
CP in jth image Ij (j = 1,2) now has covariance matrix 
{C,j/Nj^k)S where S = ®j^®). The CRLB for estimating 

the single molecule location X 2 ,f is calculated using the more 
general expression dUl. In Figure |6] the square root of its 
leading diagonal is compared to the large K standard deviation 
for the first dimension of the LRE given in iTTSI p. 6295] when 
registration is performed using the generalized least squares 
estimator. It is clear to see that the two expressions take very 
similar values, particularly for large values of K, and hence the 
close association between the CRLB expressions derived here 
and the large K results of m can be extended to the more 


general weighted covariance model. This result is general and 
not specihc to the microscopy setting. 

4) Low SNR: To demonstrate that the CRLB is an ap¬ 
propriate bound for low signal strengths we consider the 
same simulation study as in Section IVI-All but where the 
photon count associated with each control point is a uniformly 
distributed random variable on the interval [200, 700] and 300 
photons are collected for the single molecule. Figure |7] shows 
the CRLB is still appropriate in this setting. 


5) Estimating the CRLB: For the simulations studies pre¬ 
sented thus far the theoretical values of the CRLB have been 
possible to calculate due to artihcial knowledge of the true 
parameter values that form the parameter vector 9ftc- As 
this is the very thing that needs estimating the theoretical 
values of the CRLB is obviously unavailable to experimenters 
and therefore it becomes important to know how well we can 
estimate the CRLB given the estimated parameter values. 

We consider the same simulation set-up of Section IVI-A2I 
and now estimate the CRLB from (ITTI) using estimated values 
A, s, yi^F and {yi^i, ...,yi^K}, instead of true values A, s, 
X 2 ,f and ..., respectively. In Figure [8] we plot 























PREPRINT. COPYRIGHT (C) 2010 IEEE 



Fig. 6: The line indicates the square root of the CRLB for 
the first dimension of the unknown single molecule position 
parameter X 2 .f (and hence a lower bound for the standard 
deviation of the LRE) in object space dimensions, plotted 
as a function of the number of CPs. CPs are in a grid 
configuration (see Section |Vl] for more details). The crosses 
are the theoretical standard deviation of the LRE assuming the 
large K distribution given in ifTSl . 




Number of Control Points - K Number of Control Points - K 

Fig. 7: Low SNR example (see Section IVI-A4I for more 
details). The line indicates the square root of the CRLB and 
the crosses mark the sample standard deviation plotted as a 
function of the number of CPs K for the following estimated 
parameters: (a) the first dimension of the unknown single 
molecule position parameter X 2 ,f, (b) on, the leading element 
of transformation matrix A. In this simulation study CPs are 
in a grid configuration and the transformation consists of a 
rotation and translation. The vertical axis in (a) has units as 
nanometers, the vertical axis in (b) is unitless and on the scale 
X10“'*. Results are based on a sample of 10® simulations. 


the theoretical value of the CRLB for estimating X 2 ,F 7 together 
with the maximum and minimum value of the estimated CRLB 
over 10® simulations, clearly demonstrating that estimated 
parameter values can be used by experimenters to get an 
excellent estimate of the CRLB. 


B. Experimental verification 

Here we describe the experimental set up used to verify 
the theoretical results of this paper. A bead sample was 
prepared by adsorbing a dilute solution of 100-nm Tetraspeck 
microspheres (Thermo Fisher, Waltham, MA, USA) on Poly- 




Fig. 8: Left — the line indicates the square root of the true 
CRLB for the first dimension of the unknown single molecule 
position parameter X 2 .f (the leading diagonal in (fTTT i) in object 
space dimensions, plotted as a function of shear parameter A. 
The ‘x’ marks the minimum value of the estimated CRLB 
and the ‘H-’ marks the maximum value of the estimated CRLB 
taken over 10® simulations. Right — a magnification of the 
left plot for a single value of shear parameter. 


L-Lysine (Sigma-Aldrich, St. Louis, MO, USA) coated glass 
coverslip (Zeiss, Thornwood, NY, USA). A standard inverted 
microscope (Zeiss Axiovert 200) was configured with a 
63 X 1.46 numerical aperture Zeiss Plan Apochromat objective 
lens. The beads were excited by a 488nm diode laser (Toptica, 
Victor, NY, USA) and a 635nm diode laser (OptoEngine, Mid¬ 
vale, UT, USA). The emission light from the beads was split 
into two wavelength ranges, 502.5nm-537.5nm and 657.5nm- 
694.5nm, using a dichroic filter set (FF560-Di01-25x36; FFOl- 
520/35-25; FFO1-676/37-25; Semrock, Rochester, NY, USA), 
and imaged using two identical charge-coupled device (CCD) 
cameras (iXon DU897-BV; Andor, South Windsor, CT, USA). 

The imaging experiments were carried out by illuminating 
the beads with two lasers in 100ms pulse width over 599 repeat 
acquisitions. To estimate the coordinates of the beads acquired 
from each camera, we first selected region of interests (ROIs) 
containing a bead and fitted a Gaussian model using maximum 
likelihood estimation. All computations were performed using 
custom written software in MATLAB (MathWorks, Natick, 
MA, USA). 

Acquisitions 300-599 (a total of 300) were used in our 
analysis as they showed the greatest stability in photon counts 
between acquisitions and hence localization errors are consid¬ 
ered approximately iid. An example pair of images from each 
camera that need to be registered are shown in Figure |9] 

1) Verification of Gaussian distributed measurement er¬ 
rors: Throughout this paper localization errors have been 
assumed Gaussian in order to form the likelihood function 
from which the CRLBs are derived. Here, we verify this 
assumption by analysing the empirical localization estimates 
for the beads. Figure [TO] shows quantile-quantile (QQ) plots 
for the distribution of the localization estimates. The curve 
is produced by ordering the 300 independent estimates for 
either the x or y localization coordinate into increasing order 
of size. The probability of a value less than the jth ordered 
estimate (sample quantile) is approximately pj = j/301. 
The corresponding theoretical quantile of the standard normal 
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Camera 1 x-direction 


Camera 1 y-direction 


Camera 2 x-direction 


Camera 2 y-direction 




Fig. 10: QQ-plots (see Section IV-Bl for details) for the x and y coordinates of the localization estimates in both cameras. 
Each row of plots corresponds to one of two beads analysed. 



Fig. 9: Example (acquisition 599) of the pair of images of 
the bead sample to be registered, taken with the two separate 
cameras as described in Section VI-B. 


distribution is the value tj such that pj = F{tj), where F(-) 
is the cumulative density function of the standard normal 
distribution. The values fi,..., faoo are plotted on the horizontal 
axis against the ordered estimates on the vertical axes. This 
is done for x and y coordinates in both cameras for two 
separate beads (one row of QQ plots for each bead). The 
straight line indicates the ideal fit for Gaussian samples. It 
is clear the localization estimates are Gaussian distributed to a 
close approximation, except for some minor deviations at the 
distribution tails. 

2) Registration performance: Fourteen of the fluorescent 
beads that were present in the field of view for all acquisi¬ 
tions and able to be pair-matched were considered for image 
registration. One of the beads was isolated as a feature and reg¬ 
istration performed using the weighted covariance generalized 
least squares estimators for A and s mSl with 8,9,10,11,12 and 
13 of the remaining beads. Calculating the sample variance of 
the ERE in the x-direction across the 300 acquisitions, we 


compare it to the CRLB as given in ® estimated using the 
transformation parameter estimates (example values of these 
are A = and s = [1.000,1.000]'^ (each to 3 

d.p.)). Figure fTTIdisplavs the results. The four plots correspond 
to four random permutations of the beads we register with. 
It is clear that in this experiment the CRLB, to a close 
approximation, is attained. 

VII. Concluding remarks 

We have derived the CRLB for image registration per¬ 
formance under the heteroscedastic multivariate errors-in- 
variables model. Particular focus has been given to the case 
where the covariance matrices for the errors in localizing the 
CPs are all scalar multiples of a common positive definite 
matrix, a suitable model for fluorescence microscopy. Un¬ 
der this model the CRLB for estimating the location of a 
feature/single molecule has been found and is equal to the 
lower bound of the covariance matrix of the ERE, the error in 
localizing a feature/single molecule in the registered image. In 
the simplified case of that common matrix being the identity 
and the affine transformation between the pair of images being 
a scaled version of a unitary matrix, it has been shown that 
the lower bound for the covariance matrix of the ERE exactly 
matches the previously published large K expression when 
transformation parameters are estimated with the weighted 
covariance generalized least squares estimators. Therefore (fTlTi 
can now be considered to be the lower bound for the accuracy 
with which we can localize a single molecule in a registered 
image. Beyond this, it could also be used in future to develop 
strategies for the placement of the control points so that the 
estimation errors can be reduced. 

Simulations comparing the sample standard deviation of the 
transformation parameters and an element of the ERE with 
their respective theoretical lower bounds confirm that using 
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where Shh = ELi Shg = ELi 

Sff = Eti and = ELiGF^FIG,. 

Dealing with each term individually, we can write 
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Fig. 11; Experimental data results. The crosses indicate the 
square root of the sample variance of the first dimension 
of the LRE. The dots are the bootstrapped 95% confidence 
interval. The circles indicate the square root of the CRLB in 
the first dimension as calculated from ([8]l using the estimated 
transformation parameters. 


Sff+Sgg = 


0 






(14) 

where A = [Id, and Xi^k = e-d ® i = 1 ,..., d and 

k = l,...,K. 


the weighted covariance generalized least squares estimators 
for the affine parameters appears to be efficient even for low 
numbers of control points. Experimental data validates the 
theory presented. 


Appendix B 

Here we derive the EIM for the parameter vector 6ftc 
under the most general heteroscedastic model. The EIM is 
defined as 


Appendix A 

Here we derive the EIM for the parameter vector 9tc under 
the most general heteroscedastic model. We note that under the 
affine transformation assumed then S the mean vector 
for the measured feth CP locations yk = [Vi kV 1^ given 
as 

^l,k 

Axik + s 


fJ-k = 


Therefore we have 


dy-k 

dOFc 


0 Fk 
Hk Gk 


where Fk = {e^ E ® F, Gk = {eEV ® A, Hk = 
Id ® xE,Id ■ From (EJ we have 


K 

J{dTc) = ^ 

fc=l 
K 

= E 

fe=i 

Shh 

oT 
^HG 


Xk 


Hi 

Gl 



■ 0 

0 

Fk ■ 


® ^2,fc . 

. Hk 

Gk _ 


Hk ^2,kHk 
■.-1 
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Hl^EkGk 

FlnEkFk + GlnEkGk 


Shg 

Sff +Sgg 


= e{ ^G{0FTG,X2,F\yi, ■■■, VK, Uie) 

~ I OOffg 

dC{ 9 FTG, X 2 ,F\yi, ■■■, yK, J/ 1 ,f) 


^Etg 


Given (l6]l this can be expressed as 

J{9ftg) = 
where J{9tg) is as given in (O and 


0 0 
0 J{9tg) 


+ Jf{6ftg), (15) 


K 






dOpTC 


^Fd9l^G 


The zeros in the right-hand-side of (fTsT i are a consequence of 
the control point localizations being independent of the feature 
location. We have the following identities 

dyi,F 


daij 

dyiE 

ds 

dyi,F 

dxik 


= -A-^P^^^^A-\x2,f -s) = -A-^P^^^Eff, 


= -A-\ 

= 0 , 
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where is a d x d matrix of zeros except for a 1 placed 
in the (i, j)th element. It follows from these identities that 


Dp = 


Dj' = 


5/x 




del 


dfj. 


l.F 


del 


-A-\ 


'hJ. 


d? ^d\ 


Define Dtt = DlVt^ ^Dt and Dft = 


Dtf = 

JF{eFTc) = 


Dp = DIj,. It follows that 
df^i,P 1 9^1,F 




depTC 
^ A-'^n 


^dei^c 

ilpA~'^ DpT 

Dxt 

0 0 


UpT 


c{epTc) = J ^{eFTc) = 


Cpp CPT 
Ctf Bpx 


Cpc 
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CcF Bct Bcc 



Btc 
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oT 

1 ^HG 
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D 
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D^rpADipA'^DpT 0 


0 


0 
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B'J'J' Bjq 


Shh Shg 

Bct Bcc 


sic Spp + ScG 



Cjj Cjc 




Cct Ccc 



We are interested in the term Cpp whose diagonals are the 
CRLB for the localization of the feature/molecule in image 
22 - It follows that 


Cpp = 


D^n^ IDt and 


— [DpT, 0 ] 


where 


Dtt + Shh 


oT 

^HG 


Shg 

Spp+S gg 


UpT 


= A 


^A ^ — DptMiiD'pt'^ 


-1 


■ Mil 

Mi2 


Dtt + Shh 

Shg 

M 21 

M 22 


oT 

L ^HG 

Spp+S gg 


The zeros in the final row and column can be interpreted as 
arising because estimating the feature location in Ii occurs 
before registration and therefore has no dependence on the 
CP locations. The expression in (|2| follows from (fTSl) . 


Appendix C 

We consider the CRLB block matrix 


We therefore recognise that we can write Cpp as in ®. 

Appendix D 

Here we derive Ctt^ the CRLB matrix for estimating the 
transformation parameters, under Assumption I. Consider each 
term in 0 with Hi o = al g/ 2 , ^ 2,0 = cr| q /2 and A = <,R, 
where R is a unitary matrix (rotation/reflection) and c C R’*’ 
is a scaling factor. Then from (fTSli . (fOT l and (fT4li we have 


Shh = 
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E 

k=l 


1 

Ak 


Xk 

0 

Xpk 


0 

Xk 

X2,k 


xlk 

xlk 


where Xk = 
Shg = 


where we initially use the notations Btt, Btc, Bqt and 
Bcc to distinguish these from the matrices Ctt, Ctc, Cct 
and Ccc considered in (|4|l and 0. 

We note that 
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This gives 
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It is straightforward to show that Dlj.AD.i^pA^DpT = Dtt 
and hence 
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recovering the inverse FIM from Section |III] in which only the 
transformation parameters and CP locations are considered, 
an expected result stemming from the fact that the feature has 
no involvement in estimating either the parameters and CP 
locations. 


equals 
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and 


Ctt = (^Shh — Shg (Sff + Sgg) ^ Sj^G^j 


-1 


-rf 
-ri^ 


becomes 

Ctt = 0 + cri^o) 

■ ^--1 0 
X 0 

_-ri -r 2 7 - 1/2 + 7-1 (riXf+ r 2 Xj) _ 

where 7 = (l/K) Y^^=i ^ ^ “ ')~^xix{ where S = 

and xi = ( 1 ///) Ef=i % 

and Ti = where Xi = Yl,k=i Vk^^i,k = 4 *^ ® 

2 = 1 , 2 . 


a-i0 + p-^Kv^h - a-^(a-i + = (a + 

{K^)~^P)~^Q + P~^Kv‘^l 2 , where 0 = xi^fx^f- If o = 
\a+{K-^)-^P)-^ imdh = p-^Kv^ then {A-BD'-'^C)-^ = 
{aQ + 6 / 2 )-! = 6 - 1(6 + ar^)-i (a 0 * + 6 / 2 ), where 0 * = 
RQR^, with R the f /2 rotation matrix. This gives {A — 
BD-^C)-^ = PK-^v-^{ap-^ + Kf"^ + r^ + Kf^)-^Q* + 
c/ 2 , where c = {P~^Kv'^ + (a + P)~^r'^)~^. In a 

further condensing of notation we define P = { 1 /K)'y~^p. 
With Q*xi^f = 0 it can be shown that 

Cff = («-! — q;-^ (cr^ — 2cr^(a-i + P~^)~^a~'^ + 

(a-i + ,5-1)-! + cr^a~‘^{a~^ + P~^)~'^))~^ h 
= a(l — a-i(cr^(a-i(a-i + P~^)~^ — 1 )^+ 
(c,-l+^-l)-l))-l/2 


Appendix E 

Here we consider Cff, the CRLB for estimating the 
location of the feature/single molecule in the registered image, 
under Assumption 1 . Appendix |D] shows that 
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under the weighted covariance model and with Hi 0 = CTi 0 ^ 2 , 
fl 2 ,o = cr| g/ 2 , = CTi pl 2 and A = <;R where i? is a 

unitary matrix (rotation/reflection) and c € is a scaling 
factor, it follows that 


it follows that 

/ 


Cff — Cl 


Therefore the result follows from ([ 8 ]l. 

Appendix F 

Here we derive Cff under Assumptions I and II. With 
( 1 /^) Ef=i ^Xfe = and (l/K) vl^Xj^k = 0 
(doll becomes 
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with a = y^cr\ p and P = Ycfi^ + cr| q. 

Let A = a~^xi,Fxyp + P~^Kf'^I 2 , B = a~^xi^F, 
, D = a-’i + K-fp-'^, then A - BR-^C = 


C = 


''l.F 


= {a + P) 1 + 


a(^(a-i+ri)) 


= {a + P) 1 + ^ 


= « + / 3(1 + 5 


= a + 


1 
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